/*
 * ReactionSystem.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "ReactionSystem.h"

// C++ includes
#include <cmath>
#include <iomanip>
#include <set>

// Boost includes
#include <boost/foreach.hpp>

// Eigen includes
#include <Eigen/Core>
#include <Eigen/LU>
using namespace Eigen;

// GeoReact includes
#include <Utils/Algorithms.h>
#include <Utils/Assert.h>

ReactionSystem::ReactionSystem()
{}

ReactionSystem::ReactionSystem(const vector<Reaction>& reactions) :
reactions(reactions)
{
	set<string> species_set;
	
	BOOST_FOREACH(const Reaction& reaction, reactions)
		species_set.insert(reaction.GetSpecies().begin(),
		                   reaction.GetSpecies().end());
	
	species.assign(species_set.begin(), species_set.end());
}

const unsigned
ReactionSystem::GetNumReactions() const
{
	return reactions.size();
}

const unsigned
ReactionSystem::GetNumSpecies() const
{
	return species.size();
}

const vector<Reaction>&
ReactionSystem::GetReactions() const
{
	return reactions;
}

const vector<string>&
ReactionSystem::GetSpecies() const
{
	return species;
}

const Reaction&
ReactionSystem::operator[](unsigned idx) const
{
	return reactions[idx];
}

const MatrixXd
ReactionSystem::AssembleStoichiometricMatrix(const vector<string>& species) const
{
	const unsigned Nr = reactions.size();
	const unsigned Ni = species.size();
	
	MatrixXd uri = MatrixXd::Zero(Nr, Ni);
	
	for(unsigned r = 0; r < Nr; ++r) for(unsigned i = 0; i < Ni; ++i)
		uri(r, i) = reactions[r].Stoichiometry(species[i]);
	
	return uri;
}

const VectorXd
ReactionSystem::EquilibriumConstants(double T, double P) const
{
	VectorXd Kr(reactions.size());
	
	for(int r = 0; r < Kr.rows(); ++r)
		Kr[r] = reactions[r].EquilibriumConst(T, P);
	
	return Kr;
}

const VectorXd
ReactionSystem::ReactionQuotients(const VectorXd& a) const
{
	VectorXd Qr(reactions.size());
	
	for(int r = 0; r < Qr.rows(); ++r)
		Qr[r] = reactions[r].ReactionQuotient(a);
	
	return Qr;
}

const VectorXd
ReactionSystem::KineticRates(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	VectorXd Rr(reactions.size());
	
	for(int r = 0; r < Rr.rows(); ++r)
		Rr[r] = reactions[r].KineticRate(T, P, n, a);
	
	return Rr;
}

ostream&
operator<<(ostream& out, const ReactionSystem& reactions)
{
	BOOST_FOREACH(auto reaction, reactions.GetReactions())
		out << reaction << endl;
	
	return out;
}
